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NOMENCLATURE 


coefficients  in  eq  13 
2*21  em4ct2l 

parameter  defined  by  eq  1 8 

mass  specific  heat 
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phase  change  depth 

thermal  diffusivity 

o2/a, 


1  +  2  m 

temperature  penetration  depth 

— — ,  dimensionless  penetration  depth 
Pi  to,’ 

T(r'-r '•> 

■jlr,  M-r„l 

T  (r.-r0) 

integrated  temperature 
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—  X,  dimensionless  phase  change  depth 
*i 

h7  {Tm-  Tf)  (/  -  ?0>  _ 


Pi  M 


,  dimensionless 


T1 —  jT  jG(t')dt’,  dimension !ess  time 
PlCo[l  0 


6,  dimensionless  temperature 

*1 

penetration  depth 


Subscripts 


initial  value 

thawed  and  frozen  regions,  for  thaw  case 
fusion  value 
surface  value 
ambient  value 
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UNITS  OF  MEASUREMENT 

These  conversion  factors  include  all  the  significant  digits  given  in 
the  conversion  tables  in  the  ASTM  Metric  Practice  Guide  (E  380), 
which  has  been  approved  for  use  by  the  Department  of  Defense. 
Converted  values  should  be  rounded  to  have  the  same  precision 
as  the  original  (see  E  380). 
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APPLICATION  OF  THE  HEAT  BALANCE 
INTEGRAL  TO  CONDUCTION  PHASE 
CHANGE  PROBLEMS 


Virgil  J.  Lunardini 


INTRODUCTION 

Problems  of  freezing  and  thawing  arise  frequently 
in  such  diverse  applications  as  thermal  design  in 
permafrost  regions,  thermal  storage  of  latent  heat 
for  solar  systems,  and  the  heat  treatment  of  metals. 
One  is  often  interested  in  the  penetration  rate  of 
the  phase  change  interface,  the  temperature  field, 
and  the  boundary  heat  transfer  rates.  These  prob¬ 
lems  fall  into  the  category  of  conductive  heat  trans¬ 
fer  with  solidification  phase  change.  From  an  engi¬ 
neering  design  viewpoint,  exact  solutions  are  sought 
for  geometries  and  boundary  conditions  that  are 
simple  and  yet  representative  of  significant  systems. 
Unfortunately  the  mathematical  difficulties  are 
such  that  exact  solutions  to  this  class  of  problems 
are  limited  to  a  few  very  special  geometries  and 
boundary  conditions  (Lunardini  1981).  However, 
a  number  of  approximate  methods  have  been  de¬ 
veloped  that  can  yield  solutions  acceptable  for 
engineering  design.  This  report  describes  one  of 
these  approximations:  the  heat  balance  integral 
method. 

This  method,  which  has  been  used  with  good  re¬ 
sults  for  phase  change  problems,  involves  the  concept 
of  the  temperature  penetration  depth.  Consider 
the  semi-infinite  solid  shown  in  Figure  1 .  At  a  time 
f,  after  the  surface  temperature  has  jumped  to  Tv 
the  temperature  in  the  solid  will  be  disturbed  to  a 
depth  A '(f)  +  6(f).  Beyond  this  depth,  the  temper¬ 
ature  of  the  solid  remains  at  the  initial  temperature 
Tq  and  no  energy  is  transferred  beyond  this  point. 
The  penetration  distance  X  +  5  is  analogous  to  tire 
boundary  layer  thickness  in  fluid  mechanics.  The 
heat  balance  integral  method  is  similar  to  the  mo¬ 
mentum  integral  method  in  that  the  basic  equations 
are  satisfied  on  average  over  the  volume  of  thickness 
X(t)  +  6(f).  This  avoids  solving  the  partial  differ¬ 
ential  equation  at  each  point  within  the  domain  of 
interest 


Figure  I.  Temperature  penetration  depth. 


The  conduction  equation,  with  constant  thermal 
properties,  is 

3 2T  _  3 T 

S? "> 

where  a  is  thermal  diffusivity.  Now  this  equation  is 
spatially  integrated  over  the  distance  X{t)  +  6(f). 
Thus 

Xr6  a  X*e  97  a 
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The  left-hand  side  of  this  equation  is 
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Leibniz’s  rule  for  a  general  function  is 
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Then  the  heat  balance  integral  equation  is 


dt 
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dx 
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2l^l  =  o. 


(4) 


This  equation  is  valid  if  there  is  no  phase  change. 
Consider  the  case  of  phase  change  with  the  proper¬ 
ties  of  the  frozen  region  different  from  those  of  the 
thawed  region.  Using  the  procedure  outlined  above, 
there  will  then  be  two  integral  equations  as  follows: 


<$±  T  dX  PW 
dt  *  'f  dt  "  “i  L  3x 

3T,(0.tn 

SjTJ  =0 


(5) 


The  solution  of  a  general  problem  with  superheating 
or  subcooling  (the  initial  temperature  is  above  or 
below  the  fusion  value)  will  involve  two  coupled, 
nonlinear  differential  equations  for  the  parameters 
X and  6.  The  solution  will  normally  be  tedious  and 
often  requires  a  starting  solution  to  handle  the  singu¬ 
larity  at  the  origin.  However,  assume  that  the  initial 
temperature  is  Tf.  Then  the  problem  reduces  to 
only  one  differential  equation  since  the  penetration 
distance  X+8  is  now  identical  to  the  phase  change 


depth  X: 

± 1 
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3x  J=° 

X 

h  = 

(  T ,  (x,  t)  dx  . 

(9) 

The  heat  balance  integral  method  has  been  used 
extensively  for  single  phase  problems  (Goodman 
1958, 1964,  Goodman  and  Shea  1960,  Poots  1962, 
Lardnerand  Pohle  1961,  Bell  1978)  and  also  for  the 
much  more  complicated  two-phase  problems  (Lun- 
ardini  1980,  Lunardini  and  Varotta  1981).  The 
single  phase  problems  are  also  referred  to  as  non¬ 
subcooling  problems  since  the  initial  temperature  is 
identical  to  the  fusion  temperature. 


COLLOCATION  METHOD 


di  2 
dT 


dJX^i  +  T  dX 

0  dt  f  dt 


+ 
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dr2(x,t) 
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=  0 


where 

X 

Si  =  f  7",  (x,t)dx 
J0 

X+8 

h  =  7 2  (x,t)dx 


(6) 


The  usual  heat  balance  integral  equations  for  two- 
phase  problems  are  coupled  and  the  solution  can  be 
difficult.  A  slight  variation  of  the  heat  balance  inte¬ 
gral  method  can  be  used  to  find  an  explicit  functional 
relation  between  6  and  X  that  will  uncouple  the  equa¬ 
tions  and  simplify  the  solution. 

If  eq  5-7  are  added  together  the  result  will  be  the 
overall  energy  balance  for  the  volume  of  interest: 

~dt  (fViSi  +  d2c2^2  +  Pi®* 

+  (f>2c2  "P)ci )  7”f  X  - p2c2T0(X+8)  ] 


and  Ty(X,t)  =  T2(X,t)  -  Tf,  T2{X+8,t)  =  T0  have 
been  used.  The  energy  balance  at  the  phase  change 
interface  is 

dT,  (X,t)  3  T2(X,t) 

*1  3x  3x 


(7) 


3T,(0,t) 
=  ’  dx 


(10) 


The  term  (p2c2  -p,c,)  Tf  dX/dt,  in  eq  10,  is  the 
net  sensible  flux  of  enthalpy  at  the  phase  change 
interface  due  to  the  sudden  jump  in  the  specific  heats 
of  the  frozen  and  thawed  volumes.  This  term  was 
omitted  in  a  recent  study  by  Vuen  (1980),  although 
Yuen’s  derivations  implicitly  assumed  that  p2c2  *  p2c, 


2 


at  the  phase  change  interface.  The  retention  of  the 
sensible  enthalpy  term  gives  better  numerical  com¬ 
parisons  to  exact  solutions. 

Equation  7  can  be  rewritten  as  two  collocation 
equations  (see  Lunardini  1981): 


*T,(X,t) 

dx 


ar2(x,  t) 
dx 
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3 ^1)  j  dTx(X,t) 
3x2  /  3x 

ar^x.r)  bT2(x,  t) 

dx  +  *2  fa 
Z2T2{X,  t)  /  37-2(A-(  t) 


=  -p2ta 2 


For  semi-infinite  solids  the  following  temperature 
approximations  can  be  used: 

r,  =rf+tf1(x-Af)+o2(x-A')2  (13) 

T2  =  Tf  -  2  (x-X) 

+  (*-*)2  •  04) 

Equation  13,  representing  the  temperature  in  the 
region  which  has  changed  phase,  contains  two  un¬ 
known  coefficients.  One  of  these  can  be  found  from 
the  specified  boundary  condition  at  x  =  0.  Com¬ 
bining  eq  1 1  -1 4  yields 
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Figure  2.  Geometry  of  the  Neumann  problem. 
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Equations  1 7  and  1 8  reduce  to  those  of  Lunardini 
and  Varotta  (1981 )  for  the  single  phase  case  when 
0=0.  When  0*  0,  eq  1 7  and  1 8  agree  well  with 
the  exact  solution  but  are  less  accurate  than  the 
solution  of  Lunardini  and  Varotta  (1981);  the  max¬ 
imum  errors  of  about  15%  occur  at  low  Stefan 
numbers  with  high  Rvalues. 

Although  this  solution  is  for  a  step  change  in 
surface  temperature,  it  has  been  shown  that  the 
solution  is  valid  for  a  sinusoidal  surface  temperature 
if  the  step  change  temperature  is  the  average  value 
of  the  sinusoidal  temperature  over  one  half  cycle. 


NEUMANN  PROBLEM 

The  surface  temperature  of  the  volume  changes 
to  a  constant  Ts  at  the  start  of  phase  change  (see 
Fig.  2).  This  problem  has  beep,  solved  exactly  by 
Neumann  (c.  1 860)  and  approximately  by  Lunardini 
and  Varotta  (1981 )  using  the  heat  balance  integral. 
The  solution  to  eq  10-1 5  is 


X  =  2  i bs/o^r 


B\  +  ®21 


(a2j  +2Si)jj+C2i0+ 3c2iffi0+  -g  ®2i 

07) 


SPECIFIED  SURFACE  HEAT  FLUX 

The  problem  of  a  specified  surface  heat  flux  can 
also  be  solved  in  a  closed  form.  The  surface  temper¬ 
ature  will  increase  from  TQ  to  the  fusion  value  Tf 
when  melting  begins  (see  Fig.  3)  and  the  phase 
change  solution  can  then  be  obtained. 

The  surface  boundary  condition  is 

37\  (0,r) 

^~V-=G(f)-  <19) 

Equations  13  and  15  lead  to 

.  _  G  I  «21  *  ,1 


=  _G  °21  X  . 
Aj  S+o2lX 


3 


ical  computations.  The  surface  temperature  (for 
t  >  tQ)  is  given  by 


ot2 ,  S2  +  2  S  A 
0W  ~6m  =  "_2(A  +  a21  S)  ' 


(23) 


As  has  been  pointed  out  by  Goodman  (1964),  the 
solutions  here  are  valid  only  if  G(t )  is  nonotonically 
increasing  with  time  or  is  a  constant.  Pulse  type  heat 
fluxes  will  not  yield  correct  solutions. 


Figure  3.  Specified  surface  heat  flux  for  CONVECTIVE  SURFACE  HEAT  FLUX 

a  semi-infinite  medium. 


-  ga21 
°2  =  2Ar,  (5+a21*) 

The  collocation  method  allows  a  simple  relation  to 
be  derived  between  5  and  X.  By  using  equation  1 2 
this  is 

A=  f  +  Vt  +  “2’  85  (20) 

where 

B~  2*21  6m  +  a21  • 


A  problem  of  importance  is  that  of  heat  flow  from 
the  environment,  by  convection,  to  or  from  a  volume 
which  is  undergoing  phase  change.  The  situation  is 
shown  in  Figure  4  for  thawing.  This  problem  is 
physically  more  significant  than  the  Neumann  prob¬ 
lem  because  the  ambient  temperature  and  convective 
heat  transfer  are  specified  rather  than  the  surface 
temperature. 

The  surface  boundary  condition  is 

3T,  (0,r) 

-*i  ~x -  =h[T~  (M1  •  (24) 

Equations  1 3  and  1 5  now  yield 


Equation  1 0  can  now  be  solved  for  the  phase  change 
depth  S.  The  result  is 

S3  +  S2  |i  +a2i  0  +c21  Ojj 

[o  +  c2,  Oa+ j*2i 

+  jc21  dm  A(A.fl)=r*(A  +  Q21S). 


Again  using  eq  1 2 


(25) 


(21) 

There  is  no  exact  solution  of  this  problem  for  com¬ 
parison,  but  approximate  solutions  can  be  found  for 
the  single  phase  case  when  8m  =  0.  Equation  21  then 
reduces  to 

r*=  |  (5  +  S+>/H4S).  (22) 

This  is  exactly  the  equation  obtained  by  Goodman 
(1958)  with  the  usual  heat  balance  integral.  This 
solution  has  been  shown  to  be  in  good  agreement 
with  an  analog  solution  be  Kreith  and  Romie  (1955). 

Lozano  and  Reemsten  (1981)  derived  an  exact 
solution  for  the  single  phase  case.  The  solution  for 
Sym  =  0.2  was  essentially  identical  to  eq  22.  Un¬ 
fortunately  the  exact  solution  converges  so  slowly 
for  large  time  values  that  it  is  inefficient  for  numer- 


where 

,  _  2*2 1  6m  +  a21 


where 


The  energy  balance  equation,  eq  10,  can  now  be 
written  as 


dF  2  (0  +  ft21  o) 
dr  2 <p(o+  1)  +  a2,  (o  +2) 


(26) 
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Figure  4.  Surface  convection  for  a  semi¬ 
infinite  body. 
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Equation  26  can  be  written  as 

2r=  f  Qdo 
0 


FQ  =  (20  +  a2i  a)oSTm  +  (1  +  C210m)g 

+2[a(1  +C2)0m)  +  jC2,<?rn0l 

x  (/*  +  a2 , )  +  S-J- m  «2  + J  C210mg 

+2  |o(1  +€2,0^)  + j  C210m0]  (a+1) 

6(P  +  0t7 1  ) 

where  P  ~  <p  +  a2 1  o  and  g  ~  2  [P(1  -  0)  -  00] . 

There  is  no  exact  solution  of  eq  26  for  comparison 
but  it  can  be  shown  that  when  9m=  0  and  STm 
=  (°.»  -  0^)  =  0,  eq  26  can  be  solved  as 


r  =  +  o 


o-- 1  +yr+2r.  (30) 

Physically  this  is  a  single  phase  problem  with  the 
latent  heat  predominating.  Equation  30  is  the  quasi¬ 
steady  solution  (Lunardini  1981). 


The  numerical  solution  to  eq  28,  when  Bm  =  0, 
is  identical  to  the  heat  balance  integral  solution  of 
Goodman  (1958): 

12y  t  =  [(1  +  2y)  +  (2  +  7)0]  (1  +  ya( 2  *  0)) 


n/T 

-  47(7  -  1)  In 


+  vo(2  +  0 


1  +  \fy 

-1  +  7(2  +  o)+[T  +  yg(2  +  g)l  ^ 2 


+  (72  +  57)  %£  +  2  (y2  +  47  -  2)  o  -  1  -  27  (31 ) 

where  7=1+2  STm.  Equation  31  reduces  toeq  29 
when  STm  =0. 

Cho  and  Sunderland  (1981 )  presented  an  approxi¬ 
mate  method  of  solving  this  problem  for  the  single 
phase  case  {Bm  =  0).  Their  results  agree  very  well 
with  eq  31 ,  but  they  note  that  the  zero-subcooling 
solution  is  a  good  approximation  to  the  subcooling 
problem  when  6 m  *  0.  This  is  not  true,  as  can  be 
seen  from  the  graphs  presented  here.  The  subcooling 
has  a  very  significant  effect  upon  the  rate  of  phase 
change  and  may  be  ignored  only  at  the  risk  of  serious 
error. 

The  surface  temperature  is 

T-j  (0,t)-rf  cr(20  +  «2t  °) 

T*,-Tf  o(20  +  a21a)  +  2(0  +  a21o)  ' 

(32) 

The  nondimensional  surface  heat  transfer  rate  is 
(0  +  0t2 1  °) 

q*  =  - - -  .  (33) 

o(0  +  j  a2,o)  +  (0  +  a2,o) 

Equation  28  can  be  solved  by  simple,  numerical, 
quadrature.  Figures  5-14  are  plots  of  the  solution 
for  some  values  of  Stefan  number  and  0m,  with 
property  ratios  given  as  functions  of  the  volumetric 
water  content  for  soil  systems.  As  has  been  noted, 
the  heat  balance  integral  method  yields  solutions 
that  compare  quite  well  with  the  few  exact  solutions. 
Thus  the  graphs  presented  here  should  be  accurate 
for  normal  engineering  design,  especially  since  the 
soil  thermal  properties  will  normally  be  known  only 
to  within  10-20%. 

Storage  of  thermal  energy,  as  latent  heat,  is  be¬ 
coming  more  significant  as  solar  energy  becomes  more 
important.  In  general,  the  storage  of  thermal  energy 
will  play  an  increasingly  important  role  in  energy 
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Figure  1 1.  Surface  convection  for  soil,  xs  *  0.25, 

S 


Figure  12.  Surface  convection  for  soil,  x£  -0.50, 
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Figure  15.  Surface  convection,  B203  S 


Figure  16.  Surface  convection,  B203  S 


0.5/  /0.5  ,/  ,/  / 

,  .  /  //  ///  / 

B1-  /  //  /  //  / 

i  I  /  /  /  z'  / 

i  /  1  ///  /  33Li  F  -  67  KF 

~  /  //  / ,/  /  STm-0.05 

I  //  / //  /  [aa  Thaw 

4—)  ji  / //  S  C2t  0.82  1.22 

I  //  / //y  k2l  I  57  0.64 

u  !/'//  °2'  181  052 

fl  l/y/  - Framing 

\P^y  - Thawing 


P/ 

,  / 
(.5/ 

/ 

' 


33  Li  F -67KF 
STm-0» 

- Fr«tzin< 

- Thawing 


*r*  V 


Figure  1 7.  Surface  convection,  33  Lir-67  KF,  S7 


F/^ure  18.  Surface  convection,  33  LiF-67  KF,  S 
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Figure  19.  Surface  convection,  33  LIF~67  KF,  STm 


Flgure  20.  Surface  convection,  33  LiF-67  KF,  STm 


Table  1.  Thermal  properties  of  some  phase  change  materials.* 

Specific  Thermal 

Latent  heat  at  conductivity 

Fusion  hegt  of  Tf  at  Tf  Density 


Phase  change  temperature  fusion  (Btu/lbm°  FJ  ( Btu/hr  F  ft)  at  25°C 

material _ (°F) _ (Btu/lbm)  Solid  Liquid  Solid  Liquid  (lbm/ft*) 


B2°3 

842 

142 

0.41 

0.44 

0.9 

0.58 

115.5 

33  LiF-67  KF 

918 

266 

0.32 

0.39 

2.4-4.8 

2.30 

157.9 

67  NaF-33  MgF2 

1S30 

265 

0.34 

0.33 

2.4-4 .8 

2.69 

133.6 

12  NaF-59  KF-29  UF 

849 

257 

0.32 

0.38 

2.4-4.8 

2.60 

157.9 

♦ERDA  (1976) 


conservation  for  technically  advanced  countries. 
Figures  i  5-24  give  the  phase  change  depth  vs  time 
for  some  possible  phase  change  materials  with  the 
properties  listed  in  Table  1. 

With  these  graphs  the  phase  change  depth,  temper¬ 
ature,  and  heat  flux  can  be  predicted  as  a  function 
of  time.  The  computer  listing  is  given  for  the  numer¬ 
ical  quadrature  and  can  be  used  if  materials  with 
different  properties  are  considered. 


INSULATED  SEMI-INFINITE  BODY 


Figures  5-24  can  also  be  used  for  the  case  of  a 
slab  insulated  with  a  layer  of  material  when  the  insu¬ 
lation  temperature  is  T„,  as  shown  in  Figure  25.  The 
conductive  resistance  of  the  insulation  must  equal 
the  convective  resistance  of  the  air  layer.  Then 


d_ 

*l 


(34) 


The  dimensionless  phase  change  depth  is  then  given 
by 


°c 


(35) 
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ki 

Thowad 
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Figure  25.  Semi-infinite  body  with 
insulation  layer. 


The  graphs  Co.  then  be  used  by  assuming  that  the 
insulation  layer  has  no  latent  heat  and  phase  change 
starts  at  t  =  ta  when  the  temperature  of  the  insula¬ 
tion-slab  interface  reaches  Tf. 

The  single-phase  solution,  with  STm  =  0,  eq  30, 
can  be  rewritten  as 


V* 


f*?s  d2  + 


2M7Wf5(f-f0) 


-ki\d . 


(36) 


Equation  36  is  identical  to  the  quasi-steady  solution 
derived  by  Lunardini  (1981). 


CONCLUSION 

The  heat  balance  integral  method  can  be  applied 
to  conductive  heat  transfer  problems  with  phase  change 
to  obtain  good,  approximate,  solutions.  The  method 
is  particularly  useful  for  soil  systems  since  their  nature 
often  precludes  obtaining  accurate  data  on  the  soil 
thermal  properties.  Thus  the  use  of  approximate 
solutions  will  not  increase  the  uncertainty  of  the 
design  process. 

The  main  value  of  the  collocation  method  is  that 
it  provides  an  explicit  functional  relationship  between 
the  phase  change  depth  and  the  temperature  disturb¬ 
ance  depth.  This  relationship  will  usually  uncouple 
the  system  of  differential  equations  for  two-phase 
problems  and  can  lead  to  closed  form  solutions  or  to 
reduced  numerical  effort.  The  collocation  solution 
of  the  Neumann  problem  has  been  shown  to  be  quite 
accurate  with  a  worst  case  accuracy  of  less  than  1 5%. 
For  most  soil  systems  the  accuracy  is  within  5%.  The 
collocation  method  is  not  quite  as  accurate  as  the 
usual  heat  balance  integral  method  but  it  is  easier  to 
apply  to  two-phase  problems. 

Quantitative  values  have  been  obtained  for  the 
previously  unsolved  case  of  convection  at  the  surface 
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of  an  infinite  medium.  These  results  generalize  the 
widely  used  Neumann  solution  and  are  applicable 
to  the  same  physical  situations  as  the  Neumann 
problem. 

The  procedure  can  be  used  for  any  material  if  the 
appropriate  thermal  properties  are  supplied.  The  re¬ 
sults  of  this  report  apply  only  to  conductive  heat 
transfer  and  should  be  considered  as  first  approxima¬ 
tions  if  convection  occurs  within  the  melted  phase 
of  the  material. 
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APPENDIX  A.  PROGRAM  LISTING  FOR  NUMERICAL  QUADRATURE  OF  EQUATION  28 

This  appendix  includes  the  FORTRAN  program  for  the  numerical  quadrature  of  the  conduction 
phase  change  problem  for  a  semi-infinite  medium  with  a  convective  heat  flux  at  the  free  surface. 


C  FREEZING  CASE 


C  FREEZING  CASE 

C  CONDUCTION  PHASE  CHANGE 

C  CONVECTIVE  SURFACE  FLUX 

C  SEMI-INFINITE  MEDIUM 

IMPLICIT  DOUBLE  PRECISION  <A-H»0-Z> 

CALL  CONTRL<2«*COMOU7*,5) 

URI Tt< 1*782) 

782  FORMAT (1, ‘OUTPUT  WILL  APPEAR  IN  FILE  COHOUT*) 

WRITE<5»950  ) 

950  F0RMAT<1»5X»13H FREEZING  CASE) 

Z=100. 

CALL  T  NOU  <  •  WH  AT  K21  VALUE  WOULD  YOU  LIKE*, 29) 

READI1 »*)W 

CALL  T  NOU I  *  UH  AT  A21  VALUE  WOULD  YOU  LIKE*, 29) 

READ<1»*)A 

CALL  TNOU<  *WHAT  C21  VALUE  WOULD  YOU  LIKE*, 29) 

RE AO< 1 »  *  )  C 

CALL  TNOU ( * WH A T  STEFAN  VALUE  WOULD  YOU  LIKE*, 32) 

READ(1,*)ST 

CALL  TNOUI’HOW  MANY  THETA-M  VALULS  DO  YOU  HAVE*, 35) 

RE AD  < 1  ,*)I  COUNT 

DO  50  J= 1 , I COUNT 
WRITE <1,760) J 

FORMAT! 1,»INPUT  THETA M-«, 13) 

K  E  AD  1 1  ♦  *)  THETiM 
WRITE  <5,75 )THrTM,ST,W,C, A 
FORMAT! 1,/1,5X,10HTHFTA-M  =  ,F3.1/ 

/1,5X,VHSTEFAN  =  ,  Ffi  •  A 
/I  ,  5  X  ,  bHK  2 1  =  ,FG.2 
/ 1 , 5  X ,bHC 2 1  =  »  F  E  •  2 
/ 1 , 5  X , bH  A2 1  =  , F  6 . ?  ) 

WRITE  <5,100) 

FORMAT  < //I ,2X,5HSI GMA ,1 5X,3HT AU,1C  X,3HPHI 13X,6HSORT  TAU  ) 
TAU=C. 
c-ow=o. 

UO  125  K  =  0  »  2  0 
UPP=FLOAT<<) 

CALL  S I MP < CONST, THE T M, SIGMA ,POW,UPP*ST *Z,pH], TOTAL, ViCtA) 

sow=upp 

I AU=TAU*TOTAL 
S'JTAU  =  TAU**  (1  ./2.  ) 

WRITE  <5,700)SIGMA,TAU,PHI,SCTAU 
FORMAT  <1  ,?*»FA,l,10X»Fl?,5,EX«F12.,:>»fiX»F12,5) 

CONTINUE 
50  CONTINUE 

CALL  CONTRL<A»*COMOUT**5> 

CALL  EXIT 
END 

SUBROUTINE  SI MP < CON  IT *THETM,S IGM A , BOW ,UPP,ST ,2, PH] ,TOTAL*W,C,A) 
IMPLICIT  DOUBLE  PRECISION  <A-H*0-Z) 

D=0  • 

TOTAL=0. 

H:(  UPP  "E'OW  )  ZZ 
SIGMA=BOW*H*D 

CALL  FCT  <C0NST,THETM,S1 GMA,TOT,ST  »PHJ  ,  W , C , A ) 

TOTAL=TOT 
D=D* 1 • 


780 
7  5 

1  0  0 


7  00 
125 


V 


1 


13 


CALL  FCl<CGNSl*IHFT-'.»SIGPA«TOI«ST«PHl»i.»C«M 

701=101*4. 

TOI AL=T01AL*701 
0  =  £)*1. 

S  IGMA=H0W*H*D 

CALL  F C T  ( CONST  « THE1  »S  I  GH A  *T0T  »S1  »PH1  »W»C»t) 

101=101*2. 

101AL=J0TAL*101 

CALL0FCT<C0NS1«1H£1M«SIGMA«101*ST«PHI»>*C»A> 

101 AL= ( 1 01 AL*101>  *H/3. 

RC1URN 

END 

SU~-  ROU  1 1  NE  FC1  ICON  SI  »THtTM»SI  GMA«TGT«ST  *PH1  tU*C  »A) 

I MPLI C 1 1  OOJbLE  PRECISION  IA-H«0-Z> 

8E1=(2.*W*THE1M*A)/SI 

PH  I  =  BET* (SIGMA* I. >/?.♦< BE  1 *!|*(SIGKA*1 . >*»2 

1/4.*A*SIGHA*(S1GMA/2»*1»)*bEI  )**tl  •/?•  > 

HFLP=?«*PHI*<  SIoHA*1 • ) * A  *S  J  GK  A  * (S I GP  A*?  . > 
FUNC=(SI*SIGrtA**2*<  PH  I+l./i.*A*S  I  GPA))/(2.*PHIMSIG  BA 

1*1.) ♦A*SlGHA*(SIGMA*2.))*SIGHA*(l«*C»THhlH)*l./3» 

1  »C*THCIM*PHI 
P=(PHl*A*SIGMA)*2. 

R=S1*SIGMA*  12  .*PHI*A*S1  GMA>*  t  1  .  *  A  *  1  1  T*  ‘  '  L? 

T=2.«  <SIGHA*<  l.*C*THE7H  J*a./3.)«"  •.**»'**'  •'<t>»<PMI*  A* 

1  ( S I GMA* 1 •  )  ) 

y  =  csi*sigma**?*<i./3.)*c*thlth*help*?.'.  v  'f'^a.KjjiHrifO* 

l(l./3»)*C*lHCTM*PHI)*2.*<SIGMA*l<>>i  >NC*  (..  I  GF  A  1.1 

U  =  l  I  A*<SlGMA*l.>*rE7-»BFT*PHI)/<2.*PHl-('E  r*  ISIGMA*1.>>> 

D=2«*FUNC*<PH1*A*<SIGMA*1.)> 

q=k*t* r*u-o 

toi=u/p 

RETURN 

END 


*  0.*.  OOVIWMINT  PiUNTWO  ornct: 


Mo-«a/m 
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A  facsimile  catalog  card  in  Library  of  Congress  MARC 
format  is  reproduced  below. 
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